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We  consider  discrete  wavelet  transform  (DWT)  multiscale  products  for 
detection  and  estimation  of  steps.  Here  the  DWT  is  an  overcomplete  ap¬ 
proximation  to  smoothed  gradient  estimation,  with  smoothing  varied  over 
dyadic  scale,  as  developed  by  Mallat  and  Zhong.  We  show  that  the  mul¬ 
tiscale  product  approach,  as  first  proposed  by  Rosenfeld  for  edge  detec¬ 
tion,  is  a  nonlinear  whitening  transformation.  We  characterize  the  resulting 
non-Gaussian  heavy-tailed  densities.  The  results  may  be  applied  to  edge 
detection  with  a  false  alarm  constraint.  The  response  to  impulses,  steps, 
and  pulses  is  also  characterized.  A  general  closed-form  expression  for  the 
Cramer-Rao  bound  (CRB)  for  discrete  and  continuous-time  step-change 
location  estimation  in  independent  identically  distributed  non-Gaussian 
noise  is  developed,  generalizing  previous  results.  We  consider  location  es¬ 
timation  using  multiscale  products,  and  compare  results  to  the  appropriate 
CRB. 
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1.  Introduction 


Step  (or  edge)  detection  and  estimation  is  a  fundamental  problem  in  many 
areas  of  signal  and  image  processing  that  involves  basic  statistical  tradeoffs. 
A  classical  approach  to  this  problem  is  based  on  gradient  estimation,  noting 
that  the  gradient  will  be  large  when  a  noise-free  step  is  encountered.  This 
approach  has  been  heavily  exploited  in  image  processing;  e.g.,  see  Jain  [13, 
chapter  9].  Many  of  the  basic  image  edge-detection  operators  (see  Roberts, 
Prewitt,  Sobel)  reduce  in  1-D  to  a  finite  impulse  response  (FIR)  filter  with 
response  [-1,  0,  1].  More  general  extensions,  so-called  filtered  derivative 
methods,  combine  smoothing  with  gradient  estimation  to  reduce  noise  ef¬ 
fects.  These  methods  are  attractive  for  their  low-complexity  linear  imple¬ 
mentations.  In  addition,  they  tend  to  be  localized,  and  are  therefore  useful 
in  cases  such  as  images  that  may  be  characterized  by  relatively  high  signal- 
to-noise  ratio  (SNR)  and  multiple  closely  spaced  change  points. 

The  filtered  derivative  approach  is  motivated  by  application  of  linear  filters 
whose  output  is  maximized  when  a  step  is  encountered,  while  attempt¬ 
ing  to  maintain  good  location  estimation  (localization)  and  low  probability 
of  false  alarm.  A  filtered  derivative  method  that  has  received  a  lot  of  at¬ 
tention  is  the  derivative  of  Gaussian  (dG),  which  estimates  the  gradient 
after  smoothing  with  a  Gaussian  function.  The  dG  approach  can  be  de¬ 
rived  under  criteria  of  detection  and  localization  (see  Canny  [7],  Tagare 
and  deFigueiredo  [32],  and  Koplowitz  and  Greco  [15]).  The  problem  can 
also  be  formulated  in  terms  of  zero  crossings  of  the  second  derivative,  such 
as  the  Laplacian  of  Gaussian  approach,  which  is  equivalent  to  dG  in  1-D. 
Attempting  to  achieve  simultaneous  detection  and  estimation  results  in  a 
tradeoff  between  the  level  of  smoothing  and  the  variance  of  the  estimated 
step  location,  and  this  tradeoff  is  sensitive  to  the  edge  shape  and  SNR.  It 
is  therefore  of  interest  to  study  gradient  estimation  techniques  that  exploit 
multiple  levels  of  smoothing  (multiscale)  simultaneously. 

Many  alternative  techniques  are  based  on  formulating  the  problem  as  a 
change  in  distribution,  such  as  a  step  change  in  the  mean.  These  meth¬ 
ods  often  assume  a  moderate  to  large  sample  size  around  a  single  point 
of  change,  and  find  application  in  such  (typically  1-D)  problems  as  fault 
detection  or  vibration  monitoring;  e.g.,  see  Basseville  and  Nikiforov  [1] 
and  references  therein.  This  problem  may  be  formulated  as  one  of  quickest 
change  detection  in  an  on-line  setting;  a  Bayes-optimal  approach  was  for¬ 
mulated  by  Shiryaev  in  1961  [30].  An  important  approach  in  this  context  is 
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the  CUSUM  (cumulative  sum)  algorithm;  e.g.,  see  Basseville  and  Nikiforov 
[1,  chapter  2],  These  methods  may  be  extended  to  include  other  changes  in 
statistical  description,  such  as  changes  in  second-order  statistics.  We  note 
that  CUSUM  and  other  higher  complexity  approaches  may  outperform  fil¬ 
tered  derivative  methods  at  low  SNR  when  sufficient  data  are  available 
(e.g.,  see  Basseville  et  al.  [2-4]),  although  problems  maybe  encountered  for 
closely  spaced  change  points. 

It  is  well  known  that  wavelets  may  be  used  to  characterize  the  local  regular¬ 
ity  of  signals,  and  may  be  applied  for  singularity  detection  and  character¬ 
ization;  e.g.,  see  Mallat  and  Hwang  [21],  In  this  regard,  Mallat  and  Zhong 
[22]  developed  a  discrete  wavelet  transform  based  on  smoothed  gradient 
estimation,  where  the  level  of  smoothing  varies  with  (dyadic)  scale,  conve¬ 
niently  placing  the  earlier  work  of  Canny  and  others  in  the  wavelet  con¬ 
text.  This  implementation  is  based  on  splines,  and  closely  approximates 
the  derivative  of  Gaussian  estimator;  hereafter,  we  refer  to  the  implemen¬ 
tation  in  Mallat  and  Zhong  [22]  as  the  MZ-DWT  algorithm  (Matlab  M-file 
listings  are  provided  in  the  appendix  for  the  forward  and  inverse  MZ-DWT 
algorithm).  This  approach  is  nonorthogonal,  so  as  to  preserve  regularity  in¬ 
formation  at  each  point  in  time  for  each  scale.  Thus  the  MZ-DWT  is  simply 
composed  of  a  bank  of  FIR  filters  that  are  dyadic  in  scale,  but  whose  outputs 
are  not  downsampled.  Each  filter  corresponds  to  dG  with  various  levels  of 
smoothing.  This  has  been  applied  to  filtering  (denoising)  [21,34],  compres¬ 
sion  [22],  ECG  characterization  [18],  and  acoustic  shockwave  detection  [29]. 

A  central  question,  then,  is  how  to  process  and  combine  the  multiscale  gra¬ 
dient  information.  Mallat  and  Zhong  [21]  estimated  Lipschitz  exponents 
by  tracking  local  maxima  across  scales,  an  approach  that  is  rather  sensitive 
to  noise.  Li  et  al.  [18]  developed  ad  hoc  modifications  for  the  ECG  prob¬ 
lem.  We  also  note  the  earlier  work  of  Witkin  [33],  who  developed  ideas 
for  multiscale  characterization  of  signals.  Working  on  image  processing  be¬ 
fore  the  advent  of  wavelets,  Rosenfeld  and  coworkers  suggested  an  inter¬ 
esting  idea,  to  form  multiscale  point-wise  products  [26-28].  This  approach 
is  intended  to  enhance  multiscale  peaks  due  to  edges,  while  suppressing 
noise,  by  exploiting  the  multiscale  correlation  due  to  the  presence  of  the 
desired  signal  in  a  direct  (albeit  nonlinear)  way.  It  is  interesting  to  note  that 
Rosenfeld's  original  work  used  dyadic  scales,  both  for  ease  of  implemen¬ 
tation  and  because  this  seemed  to  work  as  well  as  any  other  combination 
of  scales.  The  multiscale  product  idea  was  recently  applied  for  filtering  of 
magnetic  resonance  images  by  Xu  et  al.  [34],  However,  statistical  analysis 
of  the  multiscale  product  technique  is  lacking  in  this  work. 
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In  this  report,  we  study  the  multiscale  product  approach  for  edge  detec¬ 
tion  and  estimation.  We  characterize  it  statistically  and  evaluate  perform¬ 
ance  for  detection  and  estimation  of  edges.  First,  we  briefly  introduce  the 
MZ-DWT  and  the  multiscale  product;  numerical  examples  throughout  the 
report  are  based  on  the  MZ-DWT  implementation.  Then  we  develop  the 
second-order  statistics  of  multiscale  products,  showing  the  whitening  ef¬ 
fect  of  this  nonlinear  operation.  We  consider  the  step  and  pulse  response 
as  a  function  of  step  width,  quantifying  the  suppression  of  impulses  in  the 
multiscale  product.  We  then  characterize  the  first-order  probability  density 
function  (pdf)  of  multiscale  products,  showing  that  they  are,  in  general, 
heavy-tailed.  The  resulting  pdfs  are  applied  to  detection  of  edges  with  a 
false  alarm  constraint.  We  consider  step  location  estimation,  and  give  re¬ 
sults  based  on  a  general  closed-form  Cramer-Rao  bound  (CRB)  derivation, 
provided  in  the  appendix.  We  close  with  a  discussion  of  our  study. 
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2.  A  Non-Orthogonal  DWT  and  Multiscale  Products 


In  this  section  we  provide  background  and  a  quick  overview  of  the  MZ- 
DWT  algorithm.  Matlab  M-file  listings  are  provided  in  the  appendix.  We 
then  introduce  the  multiscale  product  and  provide  motivating  examples. 

2.1  Multiscale  Gradient  Estimation 

Consider  a  continuous  time  wavelet  given  by  0s(f)  =  (1  /s)0(f/s),  such 
that  0(f)  meets  appropriate  criteria  to  be  a  wavelet.  The  continuous  wavelet 
transform  (CWT)  of  x(t)  £  L2(R)  is  Wsx(t)  =  x(t)  *  0s(f),  where  *  denotes 
convolution.  The  CWT  is  invertible  via 

roc  roc  We 

x(t)  =  C  /  Wsx(t)  0*(t  -  t)dr— ,  (1) 

■JO  J  —  oo  S 

where  C  is  a  constant  that  depends  on  0(f),  and  0*(f)  denotes  complex 
conjugation.  In  this  report  we  restrict  ourselves  to  1-D;  full  details  and  2-D 
extensions  of  the  CWT  are  described  by  Mallat  [21,22]. 

Here  we  are  concerned  with  wavelets  such  that  0(f)  =  du(t)/dt ,  with  u(t) 
acting  as  a  local  average  or  smoothing  function.  Now, 

Wsx(t)  =  x(t)  *  0)  =  S~[JX*  u0{t),  (2) 

provided  that  0(f)  is  a  valid  wavelet  with  zero  first  moment,  so  that  differ¬ 
entiation  and  integration  may  commute  [22],  Thus,  for  appropriate  choice 
of  w(f),  Wsx(t)  can  be  interpreted  as  a  derivative  of  a  local  average  of  *(f) 
where  the  degree  of  smoothing  depends  on  scale  s.  Of  particular  interest 
here  is  the  case  of  u(t)  closely  approximating  a  Gaussian  function;  the  re¬ 
sult  is  derivative  estimation  at  various  levels  of  smoothing. 

Mallat  and  Zhong  [22]  developed  a  discrete  wavelet  transform  (DWT)  based 
on  a  discrete-time  approximation  to  a  Gaussian  u(t)  using  a  cubic  spline. 
Consequently,  0(f)  is  a  quadratic  spline.  Restricting  to  dyadic  scales,  the 
DWT  of  x(n),  1  <  n  <  N,  consists  of 

W2jx(n),  j  =  1,2,---,  J—  1  ,  (3) 

where  J  =  log2  N,  plus  the  remaining  coarse  scale  information  denoted 
by  Sj(n).  This  DWT,  consisting  of  J  x  N  points,  is  overcomplete  (non- 
orthogonal).  This  contrasts  with  the  (perhaps  more  commonly  encountered) 
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orthogonal  wavelet  transforms  where  the  number  of  coefficients  decreases 
with  scale.  The  inverse  DWT  may  be  readily  computed,  enabling  filtering 
and  reconstruction  (see  remark  1  at  the  end  of  this  section). 

To  emphasize  the  linear  filtering  aspect  of  equation  (3)  and  to  simplify  no¬ 
tation,  let 

yj(n )  =  W2jx(n)  =  ^  hj(k)x(n  —  k)  (4) 

k 

denote  the  DWT  response  to  x(n)  at  the  j th  scale,  where  hj(n)  is  the  im¬ 
pulse  response  (IR)  of  the  jth  DWT  filter.  The  MZ-DWT  may  be  expressed 
in  the  Fourier  domain  as  follows.  The  Fourier  transform  of  <j)(n)  is  [22] 


=  ju 


sin  uj/ 4 
uj/ 4 
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(5) 


Denoting  the  Fourier  transform  of  hj(n)  as  Hj  ( eJuJ ) .  then 


'  G(u), 

3  =  1 

Hj(en  =  < 

G(2uj)H(uj), 

3  =  2 

(6) 

.  G{2^~1uj)H(2^-2uj)  ■  ■  ■  H(u): 

Go. 

V 

to 

where 

G{u)  =  Aje^/2  sin(w/2), 

(7) 

and 

H(uj)  =  ejuj/2  (cos(cj/2))3. 

(8) 

G(uj)  and  H( uj)  give  rise  to  a  time-domain  recursion  in  the  MZ-DWT  filter 
coefficients  [22],  Note  that  /T,  (eW)  is  purely  imaginary  and  hence  hj(ri)  is 
anti-symmetric. 

Hereafter  we  restrict  our  attention  to  Mallat  and  Zhong's  implementation 
(MZ-DWT),  although  our  results  are  general  for  a  family  of  linear  deriva¬ 
tive  estimation  filters.  Here,  each  filter  is  an  approximation  to  dG,  with 
smoothing  increasing  with  scale.  The  IR's  for  the  first  five  scales  of  the  MZ- 
DWT  are  shown  in  figure  1  (a)  through  (e),  and  some  frequency  response 
magnitudes  of  equation  (6)  are  shown  in  figure  2.  Note  that  the  linear  re¬ 
gions  around  the  origin  in  figure  2  correspond  to  the  frequency  response  of 
a  differentiation  filter. 


2.2  Multiscale  Products 

Working  before  the  advent  of  the  wavelet  framework,  Rosenfeld  and  cowork¬ 
ers  suggested  forming  multiscale  point- wise  products  [26-28].  This  is  in¬ 
tended  to  enhance  multiscale  peaks  due  to  edges,  while  suppressing  noise. 
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by  exploiting  the  multiscale  correlation  due  to  the  presence  of  the  desired 
signal.  The  multiscale  product  of  the  DWT  outputs  is  given  by 

h  h 

P(n)  =  [I  W2jx(n)  =  [I  yj(n).  (9) 

i  ja  3= jo 

As  further  motivation,  the  number  of  local  maxima  due  to  noise  decreases 
quickly  with  scale  due  to  the  increased  smoothing.  Specifically,  for  x(n ) 
white  Gaussian  noise  the  average  number  of  local  maxima  at  scale  s  =  2J+1 
is  half  that  at  scale  s  =  2J  [21],  Thus,  while  maxima  in  W2jx(n)  due  to  edges 
in  x(n)  will  tend  to  propagate  across  scales,  the  maxima  due  to  noise  will 
tend  not  to,  so  that  p(n)  will  tend  to  reinforce  the  signal  response  and  not 
the  noise. 

An  example  is  shown  in  figure  3,  depicting  the  noise-free  step-response 
behavior  of  the  MZ-DWT.  Shown  are  (a)  the  time  series,  (b)  through  (e) 
the  MZ-DWT  for  the  first  four  scales,  (f)  p(n )  for  jo  =  1/  and  (g)  ji  =  2 
and  3.  The  peak  in  the  DWT  at  various  scales  corresponds  to  the  step  edge. 
Because  the  peaks  align  across  the  first  few  scales,  the  product  p(n)  exhibits 
a  corresponding  peak. 

A  second  example  of  the  cross-scale  product  is  given  in  figure  4.  The  sig¬ 
nal  is  x(n)  =  As(n)  +  v(n),  with  A  =  10,  v(n)  unit-variance  white  Gaus¬ 
sian  noise,  and  s(n )  taking  on  values  of  1  in  the  ranges  [51, 150],  [201, 250], 
[301, 310],  and  [361,  365],  as  well  as  impulses  given  by  5{n  —  416),  5(n  —  442), 
and  S(n  —  445).  The  SNR  is  20  dB,  defined  as  10  logio(A2/cr2).  Here,  peaks 
do  not  align  across  arbitrarily  high  scales  because  neighboring  peaks  in¬ 
terfere  due  to  lengthening  filter  responses.  The  peaks  in  p(n)  are  generally 
well  pronounced,  except  for  those  corresponding  to  the  isolated  impulses 
between  n  =  400  and  450  in  the  original  time  series,  where  smoothing  leads 
to  weakened  response  at  the  higher  scales  (note  also  fig.  1  in  this  regard). 
This  suppression  of  impulses  and  the  effect  of  pulse  width  is  explored  in 
section  4.  Peaks  in  figure  4(f)  are  generally  positive  going  because  of  the 
even  number  of  products,  whereas  those  in  figure  4(g)  are  bipolar  and  pre¬ 
serve  the  edge  up /down  direction  information.  In  the  following  sections 
we  develop  properties  and  analyze  behavior  of  p{n). 

Remarks 

1.  The  inverse  MZ-DWT  may  be  readily  accomplished  [21,22],  In  ad¬ 
dition,  reconstruction  from  local  maxima  via  alternating  projections 
may  also  be  attained,  leading  to  edge-detection-based  filtering  [8]. 
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Figure  3.  Noise-free  step  i 
response  of  MZ-DWT:  (a)  0.5 
(a)  time  series,  (b)  0 

through  (e)  first  four  2 

scales  of  DWT,  (f)  P2(n),  (b)  1 
and(g)p3(n).  0 


2 

(c)  1 
0 

( 

2 

(d)  1 
0 

( 

2 

(e)  1 
0 
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2 

(f)  1 
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4 

(g)  2 
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Figure  4.  MZ-DWT  20 

example  showing  (a)  (a)  0 

time  series,  (b)  through  -20 

(e)  first  four  scales  of  20 

DWT,  (f)  normalized  ^  0 

P2(n),  and  (g) 
normalized  p3  (n) . 


-20 

20 


2.  The  region  of  support  (RoS)  of  hi+\(n)  is  two  plus  twice  the  RoS  of 
hi(n )  and  the  RoS  of  h  \  (n)  =  2.  These  are  tabulated  in  table  1;  see  also 
figure  1.  Also,  hi(n )  =  — hi(—n ),  so  that 

=  0,  for  alii.  (10) 

h 

3.  As  previously  noted,  the  problem  can  also  be  formulated  in  terms  of 
zero-crossings  of  smoothed  second  derivatives  (see,  e.g.,  Mallat  and 
Hwang  [21]).  Multiscale  products  could  then  be  searched  for  minima. 

4.  To  avoid  discontinuities  at  the  window  edges  when  computing  the 
MZ-DWT,  we  use  an  odd-symmetric  extension  of  the  data.  The  issue 
is  discussed  in  Mallat  and  Zhong  [22], 

5.  An  interesting  method  for  multiscale  estimation  is  to  use  the  iterative 
regularization  approach  [11,12],  which  might  facilitate  better  align¬ 
ment  of  multiscale  peaks  in  the  noisy  signal  case.  This  approach  is 
nonlinear,  and  analysis  along  the  lines  presented  here  is  consequently 
not  straightforward. 

6.  Denjean  and  Castanie  [9]  have  proposed  multiscale  detection  tests. 
They  form  a  multivariate  Gaussian  test  statistic  that  uses  various  scales 
of  a  continuous  wavelet  transform. 


Table  1.  RoS  of  MZ-DWT 
FIR  filters  as  function  of 
scale. 

3  14 

4  30 

5  62 

6  126 

7  254 

8  510 


Scale  RoS 
1  2 

2  6 
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3.  A  Nonlinear  Whitening  Transformation 


In  this  section  we  consider  the  first-  and  second-order  statistics  (SOS)  of 
the  DWT.  For  numerical  results  we  focus  on  the  MZ-DWT  implementation, 
while  noting  that  the  analysis  applies  more  generally  to  filter  banks  with 
no  downsampling,  as  previously  described. 

To  begin  our  study  of  the  DWT  properties,  consider  the  cross-correlation 
between  outputs  of  the  DWT  We  assume  that  x(n)  is  zero-mean  real-valued 
white  noise  with  correlation  function  E[x{n)x(n  +  m)]  =  fT^d(rn).  Using 
equation  (4),  we  obtain 

E[yi(n)yj(n)\  =  crlJ2hi(k)hj(k).  (11) 

k 

It  is  straightforward  to  obtain  the  correlation  coefficient  between  DWT  out¬ 
puts,  given  by 

.  =  E k  hj(k)hj(k) 

3  (E  fcl^M£A*^(fc2))1/2’  ' 

and  the  joint  covariance  matrix  C  of  the  DWT  outputs  at  time  n  has  (i.  j) 
element  given  by 

[C\ij  =  o2x'£hi(k)hj(k ;).  (13) 

k 

These  relations  are  readily  calculated  for  the  MZ-DWT  FIR  filters.  Corre¬ 
lation  coefficients  pij  are  shown  at  the  top  of  table  2  for  1  <  i,j  <  8,  and 
the  values  on  the  diagonal  of  C  are  shown  on  the  bottom.  Note  the  de¬ 
creasing  correlation  as  the  scales  are  separated,  with  (dyadic)  scales  a  dis¬ 
tance  three  apart  being  essentially  uncorrelated.  The  uncorrelatedness  of 
separated  MZ-DWT  outputs  follows  intuitively  from  the  plot  of  MZ-DWT 
frequency  responses  shown  in  figure  2. 

Next  we  consider  properties  of  the  multiscale  product.  Let 

Pk{u)  =  yi (n)  •  y2(n)  •  •  •  yxin)  (14) 

denote  the  if'th  product  of  the  outputs  of  the  DWT,  corresponding  to  jo  =  1 
and  j\  =  K  in  equation  (9).  For  simplicity  in  the  following  we  restrict 
ourselves  to  jo  =  1,  but  the  extension  to  the  arbitrary  case  is  straight¬ 
forward.  Assume  now  that  x(n)  is  zero  mean  and  independent  identically 
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Table  2.  Top:  Correlation 
coefficients  p,,  between 
MZ-DWT  outputs  yi(n) 
and  yj  ( n )  (scales  s  =  T 
and  s  =  2J  ),  for  white 
noise  input.  Bottom: 
main  diagonal  values  of 
covariance  matrix  C, 
with  <7%  =  1 . 


E\p4(n)\  =  74a;  crt 


Pij 


1.0000 

0.5345 

0.2097 

0.0759 

0.0270 

0.0096 

0.0034 

0.5345 

1.0000 

0.6444 

0.2791 

0.1074 

0.0395 

0.0142 

0.2097 

0.6444 

1.0000 

0.6787 

0.3019 

0.1176 

0.0435 

0.0759 

0.2791 

0.6787 

1.0000 

0.6868 

0.3075 

0.1201 

0.0270 

0.1074 

0.3019 

0.6868 

1.0000 

0.6888 

0.3089 

0.0096 

0.0395 

0.1176 

0.3075 

0.6888 

1.0000 

0.6893 

0.0034 

0.0142 

0.0435 

0.1201 

0.3089 

0.6893 

1.0000 

[C]u 

3.5556 

1.1905 

0.3236 

0.0825 

0.0208 

0.0052 

0.0013 

distributed  (iid).  The  mean  is  given  by 

E[pk{ji)]  =  E[yi(ri)---yK(n)]  (15) 

=  Y  ---Y  hK{li<)E[x(n  -  h)  ■  ■  ■  x(n  -  lK)\- 

h  Ik 

It  follows  that  E[pi(n)\  =  0,  and  E\p^{n)\  =  0  due  to  the  anti-symmetry 
of  hj(n )  (see  eq  (10));  the  latter  also  follows  if  x(n)  is  symmetrically  dis¬ 
tributed  such  as  Gaussian.  More  generally,  E\jjx{n)]  =  0  for  K  odd  if  the 
odd  order  cumulants  of  x(n )  are  zero. 

In  general,  E[px{n)]  /  0  for  K  even.  For  K  =  2,  invoking  the  whiteness 
assumption  on  x(n)  yields,  using  table  2, 

E[p2(n)\  =  crl  "Y  h\(i) /12(f)  =  1.1905cr^.  (16) 

i 

For  K  =  4,  expanding  the  fourth-moment  in  equation  (15)  yields 


Y  M0M0M0M0]  +  [n2(0)r34(0)  +  ri3(0)r24(0)  +  ri4(0)r23(0)  , 

(17) 

where 

rij(T)  =  Y  hi(l)hj(l  +  r)  (18) 

1 

is  the  deterministic  correlation  of  the  filter  IRs  and 


E[x\n)\  o 
74x  E[x2{n)}2 


(19) 


is  the  normalized  kurtosis.  Evaluation  of  equation  (17)  for  the  MZ-DWT 
yields 

E\pA(n)}  =  0.0045  74*  a*  +  0.4940  a4  =  0.0045  a 4  [7^  +  110.50],  (20) 
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Note  from  equation  (20)  that  unless  74^  is  large,  the  SOS-related  terms  dom¬ 
inate  the  mean  for  K  =  4.  If  x(n)  is  Gaussian,  then  74,,.  =  0  and  equation  (20) 
simplifies  accordingly. 

We  find  the  autocorrelation  of  p/dri)  as  follows.  We  have 


rPK(r)  =  E\pK(n)pK{n  + t)\ 

\K  j  \  K  l 

=  E  n  H  hi{li)x(n  -  li)  n  hj(mj)x(n  +r-  mj) 

i  1  \  h  J  j  1  \mi 

=  X]  ^1(^1) ' '  '^2  hf<(h<)  ^2  ' '  ‘5 2  h'K{m,K) 

h  lK  mi  mK 

K  K 

x  E  I  U  x(n  —  k)  ]^[  x(n  +  r  —  rrij ) 
i=i  j= 1 


(21) 


To  evaluate  equation  (21),  let  us  now  assume  that  x(n )  is  Gaussian,  so  that 
the  expectation  in  equation  (21)  is  a  (2/i  )th  joint  moment  of  Gaussians.  This 
even-order  moment  reduces  to  a  sum  of  1  x  3  x  •  •  •  x  (2 K  —  1)  terms,  where 
each  term  is  composed  of  appropriate  products  of  rt](0)  and  rtj  ( r ) ,  defined 
via  equation  (18).  For  K  =  2  the  expectation  in  equation  (21)  reduces  to 
three  terms,  and  using  the  whiteness  of  x(n)  we  have 

rP2{r)  =  <r*[rl2(  0)  +  rn(r)r22(r)  +  r2i(r)n2(r)].  (22) 

Similarly,  for  K  =  3,  the  expectation  reduces  to  a  sum  of  15  product  terms 
in  77,  (7-),  given  by 


rp3(r)  =  \r  12(0)  •  r-3i(r)  •  r23(0)  +  ri2(0)  •  r32(r)  •  n3(0)  +  n2(0)  •  r33(r)  •  ri2(0) 

+  n3(0)  •  r2i(r)  •  r23(0)  +  ri3(0)  •  r22(r)  •  ri3(0)  +  n3(0)  •  r23(r)  •  ri2(0) 

+  rn(r)  •  r23(0)  •  r23(0)  +  rn(r)  •  r22(r)  •  r33(r)  +  rn(r)  •  r23(r)  •  r32(r) 

+  ri2 (t)  ■  r23( 0)  •  ri3(0)  +  r12(r)  ■  r2i(r)  ■  r33(r)  +  r12(r)  ■  r23(r)  •  r3i(r) 

+  r13(r)  ■  r23( 0)  •  ri2(0)  +  ri3(r)  •  r21(r)  ■  r32(r)  +  n3(r)  •  r22(r)  •  r3i(r)].  (23) 

The  number  of  terms  grows  quickly,  making  the  algebra  somewhat  tedious; 
for  example,  K  =  4  results  in  105  terms. 

For  the  MZ-DWT,  rPK  (rn)  may  be  easily  calculated  for  the  various  cases. 
Plots  of  rp2(r)  and  rP3(r)  are  shown  in  figure  5.  Note  the  dominance  of 
rP2( 0)  and  rP3(0).  For  K  >  3,  tPk{t)/vPk{ 0)  ss  5(r).  This  effect  becomes 
more  pronounced  as  K  is  increased,  as  confirmed  by  analysis  and  simu¬ 
lations.  Despite  both  the  temporal  dependence  in  each  of  the  DWT  time 
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series  yk(n),  as  well  as  the  cross-correlation  between  these  outputs,  pR-(n) 
is  effectively  whitened  for  K  >  3.  This  is  intuitively  apparent  from  figure 
2,  because  the  time  domain  product  corresponds  to  convolution  in  the  fre¬ 
quency  domain. 

Figure  5.  Theoretical  (a)  8 - , - , - , - , - , - , - , - , - , - 

correlations  of  (a)  p2  (n) 

(with  bias  removed)  6  “ 

and  (b)  p3(n),  depicting  „  4_ 
whitening  effect  of  "T? 

multiscale  product.  2 

0<  > 

-2 1 - ' - ' - 1 - ' - ' - ' - 1 - ' - ' - 

-50  -40  -30  -20  -10  0  10  20  30  40  50 
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Remarks 

1.  In  deriving  rPK{r),  we  can  handle  the  case  of  x(n)  non-Gaussian  in 
a  manner  similar  to  that  above  by  expressing  the  moment  in  equa¬ 
tion  (21)  in  terms  of  cumulants  and  exploiting  the  assumed  iid  nature 
of  x(t) .  This  will  also  generally  involve  higher-order  analogs  of  equa¬ 
tion  (18).  For  example,  for  K  =  2  with  x(n)  iid  non-Gaussian,  we 
obtain 


rp2(r)  =  ^[^f2(0)  +  ni(r)r22(T)  +  r2i(r)ri2(r)]  +74x  ^ 


X>!(0M0M*  +  T-)M*  +  r) 

.  i 


(24) 


which  corresponds  to  the  Gaussian  case  solution  of  equation  (22) 
plus  an  additional  term  that  is  proportional  to  743;.  Evaluating  for 
the  MZ-DWT,  we  find  that  the  second  term  in  equation  (24)  is 

74x  o-4[0.7086<5(t)  +  0.3543<5(t  ±1)]. 
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4.  Step  and  Pulse  Response 


Having  shown  that  forming  multiscale  products  yields  a  whitening  trans¬ 
formation  with  respect  to  noise,  we  next  consider  the  response  to  deter¬ 
ministic  steps  and  pulses  of  varying  widths.  This  enables  us  to  charac¬ 
terize  the  SNR  in  the  multiscale  product  domain,  and  to  examine  the  ef¬ 
fects  of  smoothing  on  multiscale-product  edge  response  as  a  function  of 
the  pulsewidth.  The  latter  is  particularly  important  to  understand  before 
applying  the  multiscale  product  approach. 

It  is  natural  to  define  an  SNR  in  the  multiscale  product  domain.  Suppose 
the  signal  is  given  by  s(n)  =  Au(n  —  1)  +  x(n),  where  x(n)  is  additive  white 
Gaussian  noise  as  above,  A  >  0  a  positive  constant,  and  u(n)  is  the  unit 
step  function;  u(n)  =  0,  n  <  0,  and  u{n)  =  1,  n  >  0.  We  may  define  the 
input  SNR  in  decibels  as 

SNR,  =  101og10  (25) 

The  (noise-free)  multiscale  product  response  to  u(n  —  1)  at  the  step-change 
time  n  =  1  is 


K  K  /  \ 

PK{n)  =  yi(n)  =  AK  hi{k)u(n  -  1  -  k)  J  .  (26) 

i= 1  i=l\i  ) 


In  this  case  the  MZ-DWT  yields  pk( 0)  «  (1.33.4) for  all  scales  (the  MZ- 
DWT  filters  have  been  designed  so  as  to  yield  the  same  noise-free  step  re¬ 
sponse  amplitude  at  all  scales;  see  fig.  3).  We  can  now  define  the  output 
SNR  as 

[pk(  0)|u(n_i)]2 


SNR0  =  10  log 
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'  Pk 


(0)1 


(27) 


x(n) 


where  Pi<{0)\u(n-i)  denotes  the  peak  step  response  and 
variance  of  pain)  due  to  x(n).  Note  that 


rPK^)\x(n)  is  the 


[jW(0)L(n_i)]2 

rPK  (0)lj;(n) 


(28) 


so  SNR0  oc  (SNR,)a  and  the  Nth -order  product  nonlinearity  results  in 
amplification  of  SNR,.  The  SNR  input-output  relationship  is  plotted  in 


14 


figure  6,  for  K  =  2  and  K  =  3,  using  equations  (22)  and  (26).  The  line 
SNR,  =  SNR0  is  shown  for  reference.  There  is  SNR  gain  for  SNR,;  suffi¬ 
ciently  positive,  whereas  there  is  an  SNR  loss  for  lower  values  of  SNR,;  and 
the  SNR  gain  increases  as  SNR,  increases.  Here  (3  =  0.4014  for  I\  =  2,  and 
(3  =  0.5428  for  K  =  3,  so  the  loss-gain  thresholds  are  4  dB  ( K  =  2)  and 
1.33  dB  (K  =  3). 

The  above  describes  the  response  to  an  (essentially  infinite)  step  function. 
Next  we  consider  the  response  to  a  pulse  of  finite  width.  In  this  case  the  in¬ 
teraction  between  edges  becomes  an  issue  when  they  are  both  contained 
within  the  impulse  response  region  of  support.  Also,  the  smoothing  at 
higher  scales  will  smear  short  pulses  and  dampen  the  maxima  in  px(n),  as 
depicted  in  figure  4.  This  is  evident  from  study  of  the  MZ-DWT  IRs  in  fig¬ 
ure  1,  whose  values  in  the  neighborhood  of  occurrence  of  the  impulse  grow 
steadily  smaller  with  increasing  scale.  The  suppression  of  small  features  in 
Ph'(n)  is  a  direct  result  of  the  smoothing  at  the  higher  scales,  and  may  be 
considered  beneficial  or  harmful,  depending  on  one's  point  of  view.  These 
effects  can  be  quantified  by  examining  the  noise-free  response  to  a  pulse 
of  varying  width.  Figure  7(a)  shows  normalized  values  of  753(0),  where  the 
input  is  a  pulsewidth  of  no,  given  by  s(n)  =  u(n  —  1)  —  u(n  —  1  —  no),  for 
no  =  1,  •  •  • ,  11.  Here  £>3(0)  is  the  peak  response  to  the  leading  edge  of  the 
pulse.  We  plot  10  log10  753(0) ,  normalized  by  the  maximum  of  773 (0)  over  the 
range  of  no  -  Pulsewidths  of  no  >  6  yield  essentially  the  same  peak  response 
in  P3(n) ,  equivalent  to  the  ideal  unit-step  response.  The  case  of  no  =  1  corre¬ 
sponds  to  s(n)  =  5{n  —  1),  an  impulse.  Figure  7  indicates  a  (deterministic) 


Figure  6.  SNR 
input-output 
relationship  for  step 
input  in  white  Gaussian 
noise:  (a)  for  P2  (0) ,  and 
(b)  forp3(0).  Dashed 
line  SNR,  =  SNRD  is 
shown  for  reference. 
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Figure  7.  Edge 
suppression:  (a)  product 
of  scales  1,  2,  and  3,  and 
(b)  product  of  scales  2, 

3,  and  4. 


15-dB  suppression  of  impulses  in  pii(n),  relative  to  the  noise-free  step  re¬ 
sponse  that  is  achieved  for  a  pulsewidth  of  no  >  6.  The  figure  also  shows 
«  3  dB  suppression  of  edge  response  for  a  pulsewidth  of  2  (three  samples 
total). 

We  repeat  the  above,  now  using  the  product  of  scales  2,  3,  and  4,  with  re¬ 
sults  shown  in  figure  7(b).  In  this  case,  the  additional  smoothing  and  omis¬ 
sion  of  the  first  scale  leads  to  more  suppression  of  impulses  relative  to  the 
step  response  (about  31  dB).  However,  the  price  paid  is  that  a  pulsewidth 
of  10  or  more  samples  is  now  required  for  the  edge  response  to  match  that 
of  the  ideal  step  response,  and  pulsewidths  of  5  have  «  3  dB  suppression 
in  edge  response. 

While  the  results  of  this  and  the  previous  section  are  encouraging,  they  do 
not  tell  the  entire  story.  In  the  following  section,  we  examine  the  pdf  of  mul¬ 
tiscale  products,  and  find  them  to  be  strongly  heavy-tailed  non-Gaussian. 
This  has  implications  for  edge  detection  in  the  multiscale  product  domain. 
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5.  Multiscale  Products  are  Heavy  Tailed 


Next  we  study  the  univariate  pdf's  of  the  multiscale  products  P2(n)  and 
p‘i{n).  We  show  that  they  are  in  general  non-Gaussian  and  heavy  tailed.  We 
then  consider  detection,  and  use  the  pdf's  to  set  detection  thresholds. 

5.1  pdf's  of  Multiscale  Products 

Although  it  is  whitened  with  respect  to  its  second-order  statistics,  pA'(n) 
is  in  general  distinctly  non-Gaussian.  However,  determination  of  the  first- 
order  pdf  is  not  straightforward,  even  with  x(n)  Gaussian.  A  closed  form 
is  only  available  for  the  bivariate  case,  given  by  Miller  [24,  sect.  2.3]: 

f(z)  =  ~J\W\  e~wl2Z  K0  (\z\y/vwv£) ,  (29) 

7 r  v 

where  W  =  C~ 1  is  the  inverse  of  the  bivariate  covariance  matrix  with  ele¬ 
ments  1  <  i,  j  <  2,  and  Kq  is  the  modified  Bessel  function  of  the 

second  kind  and  order  zero.  For  the  MZ-DWT  we  find  that  for  the  first  two 
scales 


1 

- 

- 

-1 

- 

- 

3 

II 

o 

II 

3.5556 

1.1905 

— 

0.3937 

-0.3360 

1.1905 

1.3951 

-0.3360 

1.0035 

and  |  W\  =  0.2822;  so  that  in  this  case 

f(z)  =  ea336°*  /x0(0.6286  \z\).  (31) 

7 r 

For  third  or  higher  order  products,  we  must  resort  to  parametric  or  non- 
parametric  numerical  techniques,  as  we  describe  next. 

Figure  8  depicts  pdf's  for  two  cases  with  x(n)  white  Gaussian,  the  theoret¬ 
ical  pdf  for  K  =  2  from  equation  (31),  and  an  unsmoothed  histogram  esti¬ 
mate  for  K  =  3  based  on  5  x  106  samples.  Unit-variance  Cauchy  and  normal 
pdf's  are  shown  for  reference.  Note  the  skewness  in  the  K  =  2  case,  due 
to  the  positive  correlation  between  y\ (n)  and  !/•>  (n) .  and  because  K  is  even. 
In  this  case  the  positive  tail  is  heavy,  showing  a  strong  impulsive  nature 
in  p2(n)  on  the  positive  side,  with  a  much  lighter  tail  on  the  negative  side. 
Experiments  with  numerical  estimation  of  the  pdf  in  the  K  =  2  case  exhibit 
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Figure  8.  (a)  Theoretical 
first-order  pdf  of  P2  (n) 
and  (b)  estimated 
first-order  pdf  of  P3(n). 
Unit  variance  normal 
and  Cauchy  shown  for 
reference. 


excellent  agreement  with  the  theory.  In  the  K  =  3  case  the  pdf  is  symmetric 
with  tail  behavior  heavier  than  Gaussian  but  lighter  than  Cauchy. 

An  appealing  alternative  to  histogram  or  kernel-based  pdf  estimation  is  to 
employ  an  L-term  Gaussian  mixture  model,  given  by 

/(2)  =  ^2 ^exp(-^)’  (32) 

where  the  model  is  parameterized  by  S  =  [<rf ,  •  ■  • ,  af]  and 
A  =  [Ai,  •  •  ■ ,  Al_i],  with  Xl  =  1  —  The  model  encompasses 

a  broad  range  of  symmetric  zero-mean  pdf's;  e.g.,  see  Liporace  [19]. 
Approximate  maximum-likelihood  (ML)  estimates  of  the  parameters  are 
readily  obtained  using  the  iterative  expectation-maximization  (EM)  algo¬ 
rithm;  e.g.,  see  McLachlan  and  Krishnan  [23,  sect.  2.7].  As  is  well  known, 
each  iteration  of  the  EM  algorithm  achieves  nondecreasing  likelihood, 
although  local  convergence  may  be  the  result,  depending  on  initialization. 
However,  for  heavy-tailed  symmetric  pdfs  (all  mixture  terms  zero  mean, 
as  in  eq  (32))  very  good  results  can  be  obtained  for  small  sample  sizes 
with  minimal  dependence  on  initialization.  The  Gaussian  mixture  is 
readily  applied  to  detection  and  estimation  problems  in  heavy-tailed  noise 
[5,16,17], 
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Figure  9  shows  an  example  of  f(z)  for  pain).  obtained  from  1000  samples 
and  L  =  10  in  equation  (32).  The  EM  algorithm  was  iterated  100  times;  this 
was  more  than  adequate  in  our  tests.  For  initialization  we  used 

Ai  =  0.98,  A2  =  0.005,  A,:  =  0.0025  for  3  <  i  <  6,  A*  =  0.00125  for  7  <  i  <  10, 

E  =  <To[l,  5, 10, 15, 20, 100, 500, 1000,  2000, 4000], 


where  is  a  robust  estimate  of  the  variance  given  by  [20,  400] 


(33) 


where  a  (the  median  absolute  deviation)  is  the  median  of  \x(n)  —  x\,  and 
x  is  the  median  of  the  observations  x(n).  The  final  parameter  estimates  of 
f(z )  in  equation  (32)  are  shown  in  table  3.  The  overshoot  near  the  origin 
and  undershoot  in  the  estimated  pdf  tail  are  typical  of  our  experience,  ap¬ 
parently  due  to  the  sample  size.  We  have  obtained  useful  estimates  with  as 
few  as  100  samples  because  of  the  relatively  heavy  tails,  i.e.,  a  sample  size 
of  100  yields  enough  information  about  tail  behavior  to  form  meaningful 
pdf  estimates. 


Figure  9.  Estimated 
first-order  pdfs  of  P3  ( n ) 
comparing 
nonparametric  large 
sample  histogram  with 
parametric  ten-term 
Gaussian  mixture. 

(Panel  (b)  is  a  blowup  of 
panel  (a).)  Mixture 
parameters  estimated 
via  EM  algorithm:  1000 
samples  and  100 
iterations. 


19 


In  table  4  we  show  some  estimated  statistics  of  px(n)  for  1  <  K  <  4,  for  iid 
unit-variance  Gaussian  input  and  iid  unit-variance  Laplace  input,  includ¬ 
ing  normalized  skewness  and  normalized  kurtosis.  Assuming  zero-mean, 
normalized  kurtosis  is  defined  in  equation  (19),  and  the  normalized  skew¬ 
ness  is  given  by 


E[x3(n)\ 
E[x2{n)}3/2  ' 


(34) 


Table  3.  Estimated 
parameters  of  ten-term 
Gaussian  mixture 
approximation  to 
non-Gaussian  pdf  of 
P3(n)  with  white 
Gaussian  input  to 
MZ-DWT. 


l 

A  i 

of 

1 

0.0709 

0.0000 

2 

0.2429 

0.0035 

3 

0.2914 

0.1031 

4 

0.1562 

0.8872 

5 

0.0940 

3.0992 

6 

0.1107 

14.0930 

7 

0.0136 

121.0198 

8 

0.0102 

128.3395 

9 

0.0064 

128.6320 

10 

0.0036 

128.6911 

Table  4.  Estimated 
parameters  of  pk  (n)  for 
varying  K,  showing 
heavy-tailed  nature  of 
pdf's,  and  skewness  for 
K  even. 


Unit-variance  Gaussian  input 

K  =  No.  of  products 

1 

2 

3 

4 

Mean 

0.0000 

1.1899 

0.0004 

0.4948 

Variance 

3.5531 

6.3654 

10.1739 

6.9256 

Skewness 

0.0007 

2.4063 

0.0143 

10.4348 

Kurtosis 

-0.0010 

10.1278 

52.7261 

230.4220 

Unit-variance  Laplace  input 

K  =  No.  of  products 

1 

2 

3 

4 

Mean 

0.0000 

2.3805 

0.0029 

2.0214 

Variance 

2.6666 

5.8274 

10.9141 

13.3247 

Skewness 

-0.0003 

4.0264 

-0.1553 

18.0487 

Kurtosis 

1.4981 

32.0013 

171.7648 

825.0208 
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These  estimates  were  obtained  by  averaging  over  100  runs,  with  105  sam¬ 
ples  for  each  run.  They  demonstrate  the  heavy-tailed  nature  of  pK(n),  as 
well  as  the  strong  skewness  for  K  even.  These  may  also  be  evaluated  theo¬ 
retically,  although  the  task  is  tedious  as  K  increases. 


5.2  Detection  in  the  Multiscale  Product  Domain 


Knowledge  of  the  pdf  of  pK(n)  may  be  used  to  set  detection  thresholds  to 
achieve  a  constant  false  alarm  rate.  This  can  be  achieved  with  a  known  or 
estimated  noise  variance.  We  show  a  1-D  edge  detection  example  in  fig¬ 
ure  10,  where  detection  is  declared  when  |ps(n)|  exceeds  a  threshold.  The 
threshold  was  obtained  from  the  histogram  estimate  of  the  pdf  for  ps(n), 
and  set  to  achieve  Pfa  =  0.01  (the  resulting  threshold  was  13.5053).  Three 
different  edges  were  used  in  additive  white  Gaussian  noise  with  variance 
=  1 .  Consider  a  sigmoidal  step  change  signal,  given  by 


xs(n) 


m.2  +  nn\e  aT(n  n°  A 

l  -|_  e-aT(n-n0-r) 


(35) 


This  model  was  used  by  Reza  and  Doroodchi  [25].  The  parameters  m\  and 
m2  are  the  signal  levels  before  and  after  the  step,  T  is  the  sampling  interval 
time,  parameter  a  determines  the  risetime,  and  the  step  occurs  (in  contin¬ 
uous  time)  at  n^T  +  r.  Without  loss  of  generality  we  can  assume  T  =  1, 
because  T  can  always  be  incorporated  into  a.  Let  A  =  m-2  —  rn \  denote 
the  step  change.  An  alternative  expression  for  xs(n)  in  terms  of  the  tanh 
function  is 


xs(n) 


A 

2 


[1  +  tanh 


+  m  i, 


n  =  0, ...,  N  —  1  . 


(36) 


For  the  example  of  figure  10,  two  of  the  signals  were  based  on  xs(n).  The 
parameter  r  was  uniformly  distributed  in  [0,T],  modeling  the  effects  of 
time  quantization  via  sampling.  The  SNR  is  defined  in  decibels  as 

SNR=  101og10^  .  (37) 

Two  cases  of  equation  (36)  were  implemented,  a  =  4.4  and  a  =  1,  cor¬ 
responding  to  a  fast  and  slower  step  risetime,  respectively.  In  both  cases 
m i  =  0  and  T  =  1,  and  m-2  was  adjusted  to  achieve  the  desired  SNR.  The 
third  curve  in  figure  10  corresponds  to  a  unit  step,  from  rri\  =  0  to  m2, 
where  m2  is  set  to  achieve  the  desired  SNR.  One  thousand  Monte  Carlo  tri¬ 
als  were  run  for  each  SNR  value,  and  p%  (n)  was  tested  at  the  position  where 
the  edge  should  produce  a  maximum.  As  we  would  expect,  slower  risetime 
results  in  poorer  detection  performance.  We  also  note  the  moderate  to  high 
SNR  required  for  a  high  probability  of  detection. 
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6.  Estimation  of  Step-Change  Location 


Next  we  consider  estimation  of  step-change  location.  This  assumes  that  a 
change  has  occurred  in  the  observation  interval,  hence  there  is  no  detection 
problem.  In  the  appendix  we  develop  a  general-closed  form  CRB  for  step- 
change  estimation  in  additive  iid  non-Gaussian  noise,  and  provide  specific 
expressions  for  the  step-change  signal  model  given  by  equation  (35)  or  (36). 
This  generalizes  the  discrete-time  results  of  Reza  and  Doroodchi  [25]. 

The  noisy  discrete-time  signal  is  given  by 

y(n)  =  x(n)  +  v(n)  =  s(a(n  —  n0  —  r))  +  v(n),  n  =  0, ...,  N  —  1  (38) 

where  v(n)  is  zero-mean  white  Gaussian  noise  with  variance  <7^.  As  de¬ 
scribed  in  the  previous  section,  the  step-change  is  centered  at  time  nQ  +  r, 
where  r  is  a  sample-time  ,  and  a  is  a  scaling  parameter  that  determines  the 
risetime  of  the  step  (larger  a  corresponds  to  a  faster  risetime).  A  general 
expression  for  the  Fisher  information  is  given  in  the  appendix  by  equation 
(A-8)  for  t  random  and  uniformly  distributed  over  the  sampling  interval; 
an  asymptotic  expression  is  also  given  in  equation  (A-9).  These  expressions 
indicate  that  the  Fisher  information  is  proportional  to  the  risetime  parame¬ 
ter  a ,  and  inversely  proportional  to  the  noise  variance. 

To  quantify  the  general  case,  we  consider  the  step  change  signal  model  of 
equation  (35)  or  (36).  For  this  sigmoidal  step  model  we  find  the  closed- 
form  Fisher  information  to  be  equation  (A-13);  an  asymptotic  form  is  also 
provided  in  equation  (A-14).  Figure  11  plots  the  Fisher  information  of  equa¬ 
tion  (A-13)  as  a  function  of  aTN  and  d.  where  T  is  the  sampling  time,  N 
is  the  window  size,  and  d  =  n0/N  is  the  fractional  index  of  the  location 
parameter  (see  the  discussion  in  the  appendix).  Thus  d  ~  0.5  corresponds 
to  the  step  change  occurring  approximately  at  the  center  of  the  observation 
window,  and  larger  values  of  aT N  indicate  that  the  signal  variation  is  more 
completely  captured  within  the  observation  window.  From  figure  11  we  see 
that  the  Fisher  information  reaches  its  asymptotic  maximum  for  a  range  of 
d  of  approximately  0.2  <  d  <  0.8,  for  aTN  >  15.  Even  for  aTN  =  2  the 
Fisher  information  is  about  half  its  asymptotic  value.  This  indicates  that 
only  local  information  around  the  step  change  is  necessary  to  achieve  good 
estimation  of  location. 

It  is  also  interesting  to  consider  the  case  when  r  is  a  nonrandom  fixed  con¬ 
stant  tq.  The  Fisher  information  for  this  case  is  given  in  the  appendix  in 
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equation  (A-ll).  We  plot  this  case  in  figure  12,  showing  Fisher  information 
versus  parameter  aT ,  with  curves  parameterized  by  0  <  r  <  0.5.  Due  to 
symmetry  the  case  of  t  =  0.6  is  equivalent  to  that  for  t  =  0.4,  and  so  on. 
The  case  of  r  =  0  (or  r  =  1)  corresponds  to  the  step-change  being  cen¬ 
tered  exactly  at  the  sample  time,  so  that  the  Fisher  information  increases 
without  bound  as  aT  increases.  The  curves  for  r  /  0  in  figure  12  show  the 
sensitivity  to  r  in  estimating  step  location.  Depending  on  the  step  risetime 


Figure  11.  Fisher 
information  for 
step-change  location 
estimation  with 
uniformly  distributed 
sample-time  offset  r. 


Figure  12.  Fisher 
information  for 
step-change  location 
estimation  with 
sample-time  offset  r 
nonrandom.  Curves 
parameterized  by  r; 
r  =  0  implies  sample 
time  and  step-change 
aligned,  r  =  0.5  implies 
step-change  center 
exactly  between 
samples. 
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and  the  sampling  time,  it  matters  where  the  samples  fall  on  the  continu¬ 
ous  representation  of  the  step.  We  note  that  as  aT  —  oo,  the  step  risetime 
becomes  instantaneous;  in  this  limit  all  of  the  curves  decrease  to  a  limit¬ 
ing  value  (except  the  r  =  0  and  r  =  1  cases)  as  the  samples  will  bracket  the 
instantaneous  change  and  no  sub-sample  time  accuracy  can  be  determined. 

As  an  example,  we  show  in  figure  13  the  theoretical  CRB,  as  well  as  exper¬ 
imental  mean-square  error  (MSE),  for  two  edge-location  estimation  meth¬ 
ods.  The  step  was  modeled  via  equation  (35)  with  r  uniformly  distributed. 
The  two  estimation  methods  are  first,  based  on  p^,{n),  and  second,  based 
on  a  simple  gradient  estimator  with  FIR  given  by  [—1,  0, 1].  The  sigmoidal 
function  was  generated  with  m ]  =  0,  T  =  1,  and  a  =  4.4,  corresponding 
to  a  risetime  of  about  4T  (see  Reza  and  Doroodchi  [25]).  The  step  height 
m2  was  set  to  achieve  the  desired  SNR.  We  used  N  =  256  and  averaged 
over  5000  Monte  Carlo  trials  for  each  SNR  value.  For  each  trial  r  was  se¬ 
lected  uniformly  in  [—0.5,  0.5].  The  step  was  centered  at  N/ 2,  and  the  es¬ 
timate  of  step  location  was  taken  to  be  the  maximum  in  a  32-sample  ob¬ 
servation  window  centered  on  N / 2.  The  simple  gradient  estimator  has  no 
smoothing,  whereas  pz{n)  exploits  multiple  smoothing  levels.  At  low  SNR 
the  benefit  of  the  smoothing  is  evident,  as  the  DWT  method  outperforms 
the  simple  gradient  estimator  by  about  5  dB.  At  high  SNR  both  methods 
are  time-quantization  error  dominated.  The  simple  gradient  approaches 
the  minimum  variance  of  1/12  of  a  sampling  time,  which  arises  due  to  the 
uniform  distribution  of  the  sampling  phase  error.  The  DWT  method  ap¬ 
proaches  1/2,  reflecting  the  effects  of  smoothing  via  the  increased  variance 
of  the  estimate. 


Figure  13.  MSEs  and 
theoretical  CRB  for 
step-change  location 
estimation. 
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Changing  a  changes  the  effective  risetime  of  the  sigmoidal  signal.  Decreas¬ 
ing  a  creates  a  slower  rising  step.  This  results  in  a  shift  upward  of  the  CRB 
as  the  step  location  is  inherently  more  difficult  to  estimate.  Repeating  the 
experiment  in  figure  13  with  lower  a  results  in  a  shift  to  the  right  as  higher 
SNR  is  needed  to  obtain  the  same  performance,  while  the  general  behavior 
of  the  curves  is  the  same. 

Remarks 

1.  In  forming  multiscale  products,  the  individual  outputs  y,(n)  could  be 
raised  to  different  powers  before  the  product  operation.  Alternatively, 
one  could  search  for  minima  in  the  reciprocal  multiscale  product.  An 
interesting  special  case  arises  in  products  of  the  form 

m 

Xm=  II(yi)'2j-  (39) 

3= 1 

For  Yj  iid  Gaussian,  Xrn  is  o-stable,  with  a  =  2_m,  and  skewness 
parameter  f3  =  1  [6].  With  these  values  for  a,  Xrn  has  finite  mo¬ 
ments  only  for  orders  less  than  2_m,  which  decreases  rapidly  as  m  in¬ 
creases,  resulting  in  very  heavy- tailed  distributions  for  Xm.  The  case 
of  m  =  1  is  well  known,  corresponding  to  the  Pearson  V  or  Levy 
pdf.  One  could  also  use  log  Xrn .  Applying  this  result  to  MZ-DWT 
outputs  would  require  skipping  scales  in  order  to  approximate  the 
iid  assumption  (see  table  2). 
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7.  Discussion 


We  have  analyzed  the  use  of  products  for  nonlinearly  combining  multi¬ 
scale  wavelet  outputs  for  detecting  and  estimating  steps  and  edges.  The 
smoothed  gradient  DWT  developed  by  Mallat  and  Zhong  provides  a 
wavelet  framework  for  this  approach.  The  technique  exploits  the  low  com¬ 
plexity  of  the  DWT  computation,  amounting  to  a  few  FIR  filters.  Despite  its 
nonlinear  nature  the  multiscale  product  is  whitening,  although  the  result¬ 
ing  noise  pdf  is  distinctly  non-Gaussian  and  generally  heavy  tailed.  Closed- 
form  expressions  are  unavailable  for  products  of  order  greater  than  two, 
but  pdf  estimates  may  be  obtained  in  a  relatively  easy  manner.  The  heavy¬ 
tailed  nature  of  the  multiscale  products  is  problematic  from  a  detection 
viewpoint,  as  edge  detection  in  the  multiscale-product  domain  amounts 
to  detecting  impulses  in  impulsive  noise.  This  points  to  the  more  general 
idea  that  gradient  estimation  is  not  particularly  well  posed,  and  that  these 
methods  are  likely  to  be  used  for  detection  only  in  a  moderate  to  high  SNR 
regime  as,  for  example,  in  image  processing  problems.  Location  estimation 
may  be  characterized  by  comparison  with  the  general  CRB.  The  results  in¬ 
dicate  tradeoffs  possible  in  exploiting  multiple  smoothing  levels  when  the 
most  appropriate  single  smoothing  level  is  not  known  a  priori. 
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Appendix  A.  Cramer-Rao  Bounds  and  M-File  Listings 


We  derive  a  closed-form  expression  for  the  Fisher  information  matrix  (FIM) 
for  a  general  step-change  problem,  and  give  specific  results  for  a  sigmoidal 
(and  equivalent  tanh)  step-change  signal  model.  This  leads  to  a  simple  ex¬ 
pression  for  the  Cramer-Rao  bound  (CRB)  for  step  change  location  estima¬ 
tion.  This  generalizes  the  discrete-time  results  of  Reza  and  Doroodchi  [25]. 
The  continuous-time  model,  considered  by  Kakarala  and  Hero  [14],  can 
also  be  generalized  to  an  arbitrary  signal  shape  along  the  lines  below. 

A.l  General  Formulation 

The  step-change  signal  is  parameterized  by  a  location  parameter  p  = 
T(n0  +  r),  where  T  is  the  sampling  interval,  n0  is  an  integer,  and  r  G  [0, 1) 
models  sample  quantization.  Initially,  we  will  assume  that  p  is  nonrandom. 
Since  our  discrete  wavelet  transform  (DWT)-based  (and  other)  algorithms 
yield  estimates  of  n0,  not  of  p,  we  will  later  model  r  as  being  random  and 
uniformly  distributed  over  [0,1).  The  bound  obtained  by  averaging  the  FIM 
for  p  wrt  r  may  be  loosely  considered  as  a  bound  on  the  estimate  of  n0;  it  is  a 
bound  on  the  average  performance.  Note,  however,  that  since  n0  is  integer 
valued,  its  CRB  does  not  formally  exist. 

Let  a  denote  a  parameter  that  characterizes  the  step  risetime.  The  noisy 
signal  is  given  by 

y(n)  =  x(n)  +  v(n)  =  s(a(n  —  n0  —  r))  +  v(n),  n  =  0, ...,  N  —  1  ,  (A-l) 

where  s(-)  is  the  deterministic  signal,  and  v(n)  is  the  additive  noise.  Note 
that  the  observed  process  y(n )  is  nonstationary — it  has  a  time-varying  mean 
given  by  the  deterministic  signal  x(n);  its  cumulants  are  identical  to  those 
of  v(n),  so  that  if  v(n)  has  time-invariant  cumulants,  so  will  y(n).  In  the 
sequel,  we  assume  that  v(n)  is  iid,  zero-mean,  non-Gaussian,  with  finite 
variance  af,,  and  pdf  pv(v).  The  colored  noise  case  can  be  handled  along  the 
lines  of  Ghogho  and  Swami  [10],  and  Swami  [31]. 

We  will  assume  that  (1)  the  deterministic  signal,  s(-),  is  parameterized  by  a 
finite-dimensional  parameter  vector  0S,  and  the  noise  pdf  pv(v)  by  a  finite¬ 
dimensional  parameter  vector,  6V;  (2)  the  elements  of  9S  do  not  depend 
upon  those  of  0V;  and  (3)  the  noise  pdf,  pv(v),  satisfies  the  regularity  condi¬ 
tions  for  the  existence  of  the  CRB. 
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Under  these  conditions,  it  was  shown  by  Ghogho  and  Swami  [10]  that  the 
FIM  for  the  parameter  vector  9  =  [0S.  9V\  is  block  diagonal.  Consequently, 
the  achievable  accuracy  in  estimating  9S  is  the  same,  regardless  of  whether 
9V  is  known.  This  result  is  well  known  in  the  white  Gaussian  case  (less 
so  in  the  colored  Gaussian  case);  this  decoupling  between  the  signal  and 
noise  parameters  is  intuitive,  since  the  signal  parameters  are  related  to  the 
mean,  and  the  noise  parameters  to  the  cumulants.*  Hence,  we  may  assume, 
without  loss  of  generality,  that  9V  is  known.  In  this  case,  we  can  use  the 
results  of  Swami  [31]  to  write  the  (m,  n)  element  of  the  FIM  corresponding 
to  9S  as 
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y.1  dx(k)  dx(k) 
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where  Ip  is  the  Fisher  information  for  the  location  parameter  of  the  pdf 
pv(v).  Thus,  the  FIM  can  be  written  as  the  product  of  three  terms:  the  SNR, 
a  pdf-dependent  term,  and  a  signal-related  term.  In  the  class  of  symmetric 
pdfs,  7,,()  >  1,  with  equality  being  attained  in  the  Gaussian  case. 

As  we  noted  earlier,  our  DWT-based  (and  other)  algorithms  yield  integer 
estimates  of  n0;  one  could,  of  course,  augment  this  with  an  interpolation 
scheme  to  estimate  p.  Alternatively,  one  could  model  r  as  being  random 
and  take  the  expectation  of  the  FIM  in  equation  (A-2)  wrt  r,  and  thus  obtain 
an  expression  for  the  achievable  accuracy  in  estimating  p.  We  will  assume 
that  r  is  independent  of  v(n). 

Since  the  pdf  of  r  is  not  related  to  n0  or  to  the  other  parameters  of  s(-),  the 
modified  FIM  is  obtained  as 
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(A-4) 


where  the  expectation  is  with  respect  to  the  pdf  of  the  nuisance  parameter 
r.  Consequently,  the  modified  CRB  is  given  by 


*For  any  iid  process  with  finite  moments,  the  sample  estimate  of  the  mean  is  uncorre¬ 
lated  with  that  of  the  variance;  further,  the  two  are  asymptotically  jointly  normal.  The  for¬ 
mer  is  true  asymptotically  for  any  weakly  mixing  process.  Further,  if  we  write  x(t)  as  Ac(t), 
then  the  achievable  bound  on  A  ("amplitude")  is  decoupled  from  those  on  the  parameters 
describing  c(t)  ("shape"). 
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(A-5) 


var{(9s}  >  J-1, 

which  may  be  easily  calculated  for  any  specific  s(u). 

To  proceed  further,  we  need  an  explicit  model  for  x(n)  as  a  function  of  nQ 
and  the  pdf  p(r).  We  will  assume  a  uniform  pdf  for  t  (which  is,  perhaps, 
the  most  reasonable  pdf  for  r),  and  evaluate  the  FIM  for  the  general  signal 
case. 


A.2  Uniformly  Distributed  r 

We  will  focus  on  the  FIM  for  the  location  parameter  p.  Because  r  is  uni¬ 
formly  distributed  in  [0,1),  we  have 


T  _  Jvo  f1  [  dx(n) 
a2  ^ 

V  n=o 
N-l 


i  2 


dr 


dp  J 

j\-  i  r\ 

=  /  [ s'(a(n  -  n0  -  r))]2dr 

<  Jo 


rv  „  N-l  ra(n—no) 

=  — j-  ^  a  /  [s'(u)]2du 

av  n=o  J a(n—no  —  i) 


=  X]  [G{a(n  -  n0))  -  G(a(n  -  n0  -  1))], 


N—l 


(A-6) 


n= 0 


where 

G(u)  =  J  [s' {u)]2du.  (A-7) 

Thus  we  obtain  a  closed-form  expression  for  the  FIM  given  by 

J  =  a^§[G(a(N  -no-  1))  -  G(a(-n0  -  1))].  (A-8) 

Finally,  we  may  evaluate  the  CRB  via  equation  (A-5).  Note  that  the  FIM  is 
proportional  to  the  risetime  parameter  a  and  inversely  proportional  to  the 
noise  variance,  as  may  be  expected. 

For  any  reasonable  step-approximator,  s(u),  its  derivative  will  vanish  out¬ 
side  a  more  or  less  small  interval  centered  around  no,  and  we  expect  that 
G(u)  will  attain  its  asymptotic  values  rather  quickly.  The  asymptotic  value 
is  obtained  by  assuming  that  N  is  large  and  that  no  is  not  too  close  to  either 
end  of  the  observation  window, 

Joe  =  a^r[G{+ oo)  -  G(-oo)]  .  (A-9) 

a- 
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One  could  consider  the  signal  s(t)  as  a  smoothing  filter  for  the  underly¬ 
ing  ideal  step  function;  one  could  then  try  to  find  the  s(t)  that  maximizes 
the  FIM  given  by  subject  to  a  sensible  constraint  such  as  unit  energy 
It  would  be  interesting  to  consider  the  joint  estimation  of  all  the  signal 
parameters  (a,  p,  and  A). 

In  the  case  of  Gaussian  noise,  we  note  that  equation  (A-8)  is  also  the  FIM 
for  the  estimate  of  the  nonrandom  parameter  p,  when  the  observations  are 
in  continuous  time,  i.e.,  y(t)  =  s(a(t  —  p))  +  v(t),  0  <  t  <  (N  —  1  )T,  where 
v(t)  is  band-limited  white  Gaussian  noise  with  spectral  density  rr2.  This  can 
be  established  by  mimicking  the  development  in  Kakarala  and  Hero  [14]. 


A.3  Sigmoidal  Edge  Model 


Consider  the  sigmoidal  step  change  signal  described  in  equation  (35), 

m2  +  mie_aT(n_n°_'r) 


XgA)  j  +  g-aT(n-no-r) 

Let  A  =  m-2  —  m ,  denote  the  step  change. 


(A-10) 


For  the  nonrandom  parameter  p  we  obtain 

(  T\2  A2 

J  =  a  ir.7t’°^  cosh^4  (' aT(k  -  n0  -  t )/2) ,  (A-ll) 

16  av  to 

whose  value  depends  strongly  upon  r,  but  not  so  strongly  upon  na  (unless 
n0  is  close  to  0  or  N  —  1). 


We  will  now  randomize  over  r.  Using  equation  (A-10)  we  obtain 


s(u) 
An) 
G(u ) 


A  m2  +  rri\e 


1  +  e~u 

=  Acosh_2(rt/2)/4 


A 2  1 

=  —  [tanh(u/2)  —  -  tanh3(u/2)], 
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(A-12) 


which  via  equation  (A-8)  leads  to  a  closed-form  expression  for  the  FIM 
(which  was  expressed  only  as  a  summation  in  eq  (10)  of  Reza  and  Doroodchi 
[25]), 
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where  S  =  (m2  —  mi)2 /a2  =  A2 /a2  is  the  SNR. 

The  tanh(u)  function  varies  from  —1  for  u  «  0  through  0  at  u 
for  u  »  1.  Thus,  the  value  of  the  AFIM  for  xs  (n)  is 

j  _  aTS^vo  4  _  aTS'jvo 
00  “  8  3  ““  6  ' 


0  to  +1 

(A-14) 


The  AFIM  is  proportional  to  the  SNR  and  to  the  risetime  parameter  aT. 
The  AFIM  approximates  the  FIM  extremely  well  when  NaT  >  >  1  (usually 
a  few  samples  for  moderate  values  of  aT). 

By  re-indexing  the  observation  window  from  0, iV  —  1  to  1, iV  we  can 
write  the  alternative  expression 

J  =  {aT\Slv°[f(aTN(  1  -  d)/2)  -  f(-aTd/ 2)],  (A-15) 


where  d  =  n0/N  is  the  fractional  index  of  the  location  parameter,  and 

f(u)  =  tanh(fx)  —  tanh3(ix)/3. 


A.4  M-File  Listings 

In  this  section,  we  provide  listings  of  Matlab  M-files  for  the  forward  and  in¬ 
verse  Mallat-Zhong  discrete  wavelet  transform  (MZ-DWT)  algorithm  [22, 
app  B], 

%  [WT,S]  =  fwt ( J, x)  Forward  wavelet  transform 

o, 

o 

%  J  =  #  of  scales 
%  x  =  data  vector 

%  WT  =  matrix  with  wavelet  transform 
%  S  =  coarse  low-pass  time  series  remaining 

function  [WT,S]  =  mal_fwt ( J,  x) 


N  =  length (x) ; 

M  =  2  *N; 

%  normalization  coefficients 


lambda  =  [1.5,  1.12,  1.03,  1.01]  ; 
if  J  >  4 

lambda  =  [ lambda,  ones ( 1 ,  J-4 )]  ; 

end 
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%  filter  coefficients 

H  =  [0.125,  0.375,  0.375,  0.125]; 

G  =  -1*  [-2.0,  2.0]; 

%  convolution  offsets 

Gn  =  2; 

for  j  =  1 : J-l 

z  num  =  2  ~  j  -  1 ; 

Gn  =  [ Gn ,  ( ( z  num+ 1 ) / 2 ) + 1 ] ; 
end 

Hn  =  3; 
for  j  =  1 : J-l 
z  num  =  2 “ j  -  1 ; 

Hn  =  [Hn,  ( ( znum+1 ) /2 ) +znum+2 ] ; 
end 

%  compute  the  WT  at  each  scale 

%  signal  is  odd-symmetric  periodically  extended  for  borders 

S  =  [f liplr (x) , x, f liplr (x)  ]  ; 

WT  =  []; 
figure  ( 1 ) 
for  j  =  0 : J-l 


znum  =  2~j  -  1;  %  #  of  zeros 

Gz  =  in_zeros (G, znum) ;  %  insert  zeros  into  G 
Hz  =  in_zeros (H, znum) ;  %  insert  zeros  into  H 

%  compute  wavelet  transform  at  scale  j  and  store  in  WT 

Wf  =  ( 1/lambda ( j+1 )) *conv (S, Gz ) ;  %  size (Wf ) 

Wf  =  Wf (N+Gn( j  +  1)  : 2*N+Gn ( j  +  1 ) -1 )  ; 

WT  =  [WT,  Wf ' ] ; 


%  compute  next  time  series 
S2  =  conv (S, Hz ) ; 

S2  =  S2 (N+Hn ( j+1) :2*N+Hn ( j+1) -1) ; 
S  =  [ f liplr ( S  2 ), S2 , f liplr  ( S  2 ) ] ; 
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end 


S  =  S (N+l : 2*N) ; 
return 


%  x  =  iwt (WT, S)  Inverse  wavelet  transform 

o, 

o 

%  WT  =  matrix  with  wavelet  transform 
%  S  =  leftover  low-pass  time  series 
%  x  =  reconstructed  time  series 

o, 

o 

function  x  =  mal_iwt (WT, S ) 

[N, J]  =  size(WT);  %  J=#  of  scales,  N=data  length 

%  normalization  coefficients 

lambda  =  [1.5,  1.12,  1.03,  1.01]; 
if  J  >  4 

lambda  =  [ lambda, ones ( 1 ,  J-4 )]  ; 

end 

%  filter  coefficients 

H  =  [0.125,  0.375,  0.375,  0.125]; 

K  =  [0.0078125,  0.054685,  0.171875]; 

K  =  [K,  -l*f liplr (K)  ]  ; 

%  convolution  offsets 

Kn  =  3; 

for  j  =  1 : J-l 

znum  =  2 ~ j  -  1; 

Kn  =  [Kn, ( ( znum+1 ) /2 ) +2 * znum+3 ] ; 
end 

Hn  =  2; 
for  j  =  1 : J-l 
znum  =  2 ~ j  -  1; 

Hn  =  [Hn, (( znum+1 ) /2 ) +znum+2 ] ; 
end 
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%  recursively  compute  the  inverse  WT,  proceeding  down  in 
%  scales,  signal  is  odd-symmetric  periodically  extended 

S=S  ( : )  ; 

SI  =  [fliplr (S) , S, fliplr  (S) ] ;  S1=S1(:); 
for  j  =  J : -1 : 1 

znum  =  2'(j-l)  -  1;  %  #  of  zeros 

Kz  =  in_zeros (K, znum) ;  %  insert  zeros  into  K 

Hz  =  in_zeros (H, znum) ;  %  insert  zeros  into  H 

WT  j  =  WT ( : , j ) ;  WTj=WTj(:); 

WT_ext  =  [ fliplr (WTj ), WTj , fliplr (WTj )] ;  WT_ext=WT_ext ( : ) 
A1  =  lambda ( j ) *conv (Kz , WT_ext ) ; 

A1  =  A1 (N+Kn ( j )  : 2  *N+Kn ( j  )  -1 )  ; 

A2  =  conv(Hz,Sl); 

A2  =  A2 (N+Hn ( j)  :2*N+Hn ( j) -1)  ; 

SI  =  A1  +  A2  ; 

SI  =  [fliplr (SI) SI'  , fliplr (SI)  ']  ;  S1  =  S1 ( : ) ; 
end  %  end  IWT  loop 
x  =  SI  (N+l : 2*N) ' ; 
return 


%  y  =  ins_zeros (x, n) 

o, 

o 

%  Insert  n  zeros  between  elements  of  vector  x 
%  and  return  in  y.  Used  in  discrete  wavelet  transform. 

o, 

o 

function  y  =  in_zeros (x, n) 

if  n  ==  0 
y  =  x; 
return 

end 
if  n>0 

newlen  =  (n  +  1 ) *length  (x) ;  %  length  of  y 


y  =  zeros (1, newlen) ;  %  new  filter  vector 
index  =  1 : n+1 : newlen-n;  %  indices  of  data 
y (index)  =  x;  %  insert  data  into  y 

end 

return 
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